The invention relates to an apparatus for fast imaging of eg the human body using ultrasound.
The image is made using a multi-element probe in which all or selected ones of the elements are used as transmitters. The reflected signal from the object is then measured by all elements and ultrasound beams are simultaneously focused throughout the imaging region. The image is updated with the new information every time a new element or a group of elements are used as a transmitting aperture. Hereby the image is continuously updated and can be used for probing moving structures and enhance the imaging of blood velocities.
The technological development has constantly led to improvements in ultrasound beam formers as described in [1]. The introduction of digital beam formers made the dynamic receive focusing possible. Unfortunately it is not possible to dynamically control the transmitted acoustic energy, and only a fixed focus is used in transmit.
The basic concept of focusing is to constructively add energy in the combined transmit-receive signal. Ideally, by means of the receive processing it is possible to compensate for the phase changes in the transmit. This corresponds to a dynamically focused transmit and receive imaging system. A composite image (obtained by multiple emissions, focused at different depths in transmit) is the ideal reference image.
One way to obtain a dynamic transmit focus is to use synthetic array imaging. There are three variations of the synthetic aperture imaging. 1: using a synthetic transmit aperture, 2: using a synthetic receive aperture, and 3: using a combination of synthetic transmit-receive aperture. All of these combinations have been studied, and have shown different advantages and drawbacks. The time necessary to acquire a single image Tacquire image is proportional to the number of emissions Nemission, the time necessary to record the reflected ultrasound wave from a single emission Trecord emission, the number of scan lines N1, and it is inversely proportional to the number of the parallel receive beam formers Nparallel:
Taquireimage=Tacquirescanlinexc2x7Nemissionxc2x7N1/Nparallelxe2x80x83xe2x80x83(1)
The acquisition time for an image, with a typical depth of 15 cm assuming that the speed of sound is 1540 m/s, is 200 xcexcs. If Nemission=64 and N1=Nparallel then Tacquire image=12.8 ms, which results in 78 frames/sec. For every new frame in the imaging process, the data acquired for the previous frame is discarded.
The recursive imaging method and apparatus according to the invention uses the beam formed lines from the previous frame to create a new frame after every emission. This results in Nemission=1, Tacquire image=200 xcexcs, and a frame rate of 5000 frames/sec. The invention uses a synthetic transmit aperture and preferably receives with the full aperture. The prior art synthetic transmit aperture focusing [2], [3] is presented below. Based on it, the new recursive ultrasound imaging technique according to the invention is mathematically derived and discussed below.
One of the problems in emitting with only one element is the signal-to-noise ratio, since there are physical limitations to what power can be sent with a single transducer array element. The problem was previously studied [2], [4], and a solution using multiple elements in transmit was suggested in [2]. This problem is discussed below.
The other problem is the presence of motion artifacts due to the time of acquisition for a single image. This is avoided with the invention by decreasing the number of emissions necessary to create one frame. This can be done with a sparse synthetic aperture, and the results of the decreased number of emissions is presented.
Prior Art Synthetic Array Imaging
Phased linear arrays are used for creating sector B-mode images as shown on FIG. 2. The image consists of a number of scan lines N1. The scan lines l=(1 . . . N1) have a common origin. Each of them has a different angle xcex81 with the normal vector to the transducer""s surface. For simplicity all considerations are made only in the z-x plane, and it is assumed that a single transducer element emits a cylindrical wave propagating at a constant speed c in a linear medium, as shown in FIG. 3.
The measurement situation is shown in FIG. 4. Element i with center coordinates (xi,zi) emits a spherical wave. The wave front reaches point P(xP,zP) after time                                           t            iP                    =                                                                                          (                                                                  x                        P                                            -                                              x                        i                                                              )                                    2                                +                                                      (                                                                  z                        P                                            -                                              z                        i                                                              )                                    2                                                      c                          ,                            (        2        )            
where c is the speed of sound. If all the elements are excited, they will form a pressure field in front of the transducer, which is a sum of the fields created by each of the elements. The emitted waves reach point P at different time instances depending on the element positions.
In order to align the wave fronts and thereby focus the acoustic energy in point P, the signals from the individual array elements must be appropriately delayed. In the calculation of the delays, one point from the transducer is selected as a reference point. All the delays are calculated relative to the time necessary for a sound wave to travel the distance between the reference point and the focal point. For a sector image, the center of the transducer array is the usual choice. The reference point is the center of element C(xc,zc) for the case depicted on FIG. 4. The delay for element i for focusing the energy at point P(xp,yp) is calculated by:                                           d            i                    =                                    t              CP                        -                          t              iP                                      ⁢                  
                ⁢                              d            i                    =                                                                                                                (                                                                        x                          C                                                -                                                  x                          P                                                                    )                                        2                                    +                                                            (                                                                        z                          C                                                -                                                  z                          P                                                                    )                                        2                                                              -                                                                                          (                                                                        x                          i                                                -                                                  x                          P                                                                    )                                        2                                    +                                                            (                                                                        z                          i                                                -                                                  z                          P                                                                    )                                        2                                                                        c                                              (        3        )            
FIG. 5 shows the geometry used for the simulations. The origin of the coordinate system O(0,0) lies in the middle of the physical linear array transducer. All scan-lines (l=1 . . . N1) start from the origin O and have angle xcex81 with the z-axis. The focal points lie on these scan lines and have coordinates:
xe2x80x83xf1=rfxc2x7sin"THgr"1
zf1=rfxc2x7cos"THgr"1xe2x80x83xe2x80x83(4)
where rf is the axial distance to the focal point. Dynamic focusing is obtained if rf changes in time as
rf=ck/fsxe2x80x83xe2x80x83(5)
where k is the sample number and fs is the sampling frequency. In conventional imaging systems the transmission uses a fixed rf, whereas the receive can be dynamically focused along the current scan line.
One way to create a synthetic aperture image is to emit with all the elements from the transducer array one at a time (see FIG. 6). During reception data is recorded and beam formed. After a number of emissions Nxmt all the beam formed RF lines from the separate emissions are summed to create the scan-lines that are envelope-detected and displayed.
The beam formed signal from emission n with element i for line l is sli(n)(t). The number for the emission n (0 less than n less than ∞) is relative to the beginning of the continuous imaging process. This number is relevant only for the recursive imaging and in this section a simplified notation is used: sli(t)=sli(n)(t). The time t is relative time from the emission of the pulse. The RF signal that is envelope-detected and displayed is:                                           s            l                    ⁢                      xe2x80x83                    ⁢                      (            t            )                          =                              ∑                          i              =              1                                      N              xmt                                ⁢                      xe2x80x83                    ⁢                                    s              li                        ⁢                          xe2x80x83                        ⁢                          (              t              )                                                          (        6        )                                                      s            li                    ⁢                      xe2x80x83                    ⁢                      (            t            )                          =                              ∑                          j              =              1                                      N              rcv                                ⁢                      xe2x80x83                    ⁢                                                    a                lij                            ·                              s                ij                                      ⁢                          xe2x80x83                        ⁢                          (                              t                -                                  d                  lij                                            )                                                          (        7        )            
where Nxmt is the number of transmit elements, Nrev is the number of receive elements, sij(t) is the recorded data, alij is a weighting coefficient (apodization), and dlij is the delay for image line l, when emitting with element i and receiving with element j.
The delays dlij is the sum of the delays for elements i and j calculated with formula (3) for the points in line l
dlij=dli+dlj
For the case in FIG. 5, the delays are calculated by:                               d          lij                =                                                                                                  x                    fl                    2                                    +                                      z                    fl                    2                                                              -                                                                                          (                                                                        x                          t                                                -                                                  x                          fl                                                                    )                                        2                                    +                                      z                    fl                    2                                                                        c                    +                                                                                                                x                      fl                      2                                        +                                          z                      fl                      2                                                                      -                                                                                                    (                                                                              x                            j                                                    -                                                      x                            fl                                                                          )                                            2                                        +                                          z                      fl                      2                                                                                  c                        .                                              (        8        )            
This is the basic formula for a geometrical quadrature focusing. In Eq. (8) (xf1, zf1) are the coordinates of the focal point, xj is the coordinate of the center of the j""th element and xi is the coordinate of the center of the i""th element. In transmit, the delay for element i is equal to zero, and in receive the delay for element j is equal to dlij.
The synthetic array imaging assumes that the tissue under investigation is stationary, which is often not the case in medical ultrasound. The images will, thus, be blurred and the synthetic array imaging cannot be used for velocity estimation.
It is the object of the invention to overcome this deficiency and disadvantage of the known kinds of apparatus. With the invention this object is achieved by an apparatus that continually updates the image using the current pulse emission thereby obtaining a continuous image for tracking moving structures.
The concept of the invention is as follows: A new frame n is created, based on the previous frame nxe2x88x921 and the data acquired at the last emission. The calculation procedure is derived from a synthetic array aperture imaging technique as presented and further developed.
Recursive Imaging with the Invention
A simple example for synthetic imaging is shown in FIG. 6. Here n is the number for the current emission and Nxmt is the number of transmissions between two emissions with the same element. slk(n) is the beam formed signal at emission n with element k for line l. The relation between n and k is given by:
k(n)=(n mod Nxmt)+1xe2x80x83xe2x80x83(9)
The first image can be obtained after Nxmt emissions. Equation (6) can be rewritten as:                               S          t                      (            n            )                          =                              ∑                          m              =                              (                                  n                  -                                      N                    xmt                                    +                  1                                )                                      n                    ⁢                      xe2x80x83                    ⁢                                                    a                                  k                  ⁢                                      xe2x80x83                                    ⁢                                      (                    m                    )                                                              ·                              s                                  lk                  ⁢                                      xe2x80x83                                    ⁢                                      (                    m                    )                                                                    (                  m                  )                                                      ⁢                          xe2x80x83                        ⁢                          (              t              )                                                          (        10        )            
Equation (10) shows that a new frame can be formed after any emission nxe2x89xa7Nxmt by summing the beam formed lines s(m)lym),nxe2x88x92Nxmt less than mxe2x89xa6n. For two consecutive emissions the expression is:                                                         S              l                              (                                  n                  -                  1                                )                                      ⁢                          xe2x80x83                        ⁢                          (              t              )                                =                                    ∑                              m                =                                  (                                      n                    -                                          N                      xmt                                                        )                                            n                        ⁢                          xe2x80x83                        ⁢                                                            a                                      k                    ⁢                                          xe2x80x83                                        ⁢                                          (                      m                      )                                                                      ·                                  s                                      lk                    ⁢                                          xe2x80x83                                        ⁢                                          (                      m                      )                                                                            (                    m                    )                                                              ⁢                              xe2x80x83                            ⁢                              (                t                )                                                    ⁢                  
                ⁢                                            S              l                              (                n                )                                      ⁢                          xe2x80x83                        ⁢                          (              t              )                                =                                    ∑                              m                =                                  (                                      n                    -                                          N                      xmt                                        +                    1                                    )                                            n                        ⁢                          xe2x80x83                        ⁢                                                            a                                      k                    ⁢                                          xe2x80x83                                        ⁢                                          (                      m                      )                                                                      ·                                  s                                      lk                    ⁢                                          xe2x80x83                                        ⁢                                          (                      m                      )                                                                            (                    m                    )                                                              ⁢                              xe2x80x83                            ⁢                              (                t                )                                                                        (        11        )            xe2x80x83Sl(n)xe2x88x92Sl(nxe2x88x921)=ak(n)xc2x7slk(n)(n)(t)xe2x88x92ak(n)xc2x7slk(nxe2x88x92Nxmt)(nxe2x88x92Nxmt)(t)
Sl(n)=Sl(nxe2x88x921)+ak(n)xc2x7slk(n)(n)(t)xe2x88x92ak(nxe2x88x92Nxmt)xc2x7slk(nxe2x88x92Nxmt)(nxe2x88x92Nxmt)(t)xe2x80x83xe2x80x83(12)
From equation (9) it can be seen that k(n)=k(nxe2x88x92Nxmt). In recursive imaging a new frame is created at emission n by adding the new information Info(n) to the image Image(nxe2x88x921) and subtracting the information obtained at emission nxe2x88x92Nxmt, Info(nxe2x88x92Nxmt). The number of summations per sample is decreased from Nxmt in Eq. (6) to only two.
Since a new frame consisting of a number of simultaneously beam formed lines is created at every pulse emission, the image is updated after every 200 xcexcs when c=1540 m/s, and the data is acquired to a depth of 15 cm.
Add-only Recursive Imaging
To decrease the amount of the necessary storage memory the calculation procedure can be modified to the add-only recursive imaging. Consider the following equations:                     B        =                              ∑                          i              =              1                        N                    ⁢                      xe2x80x83                    ⁢                      A            i                                              (        13        )            xe2x80x83Al=A=const.
                    B        =                              N            ·            A                    =                                                    N                ·                A                            +              A              -              A                        =                                                                                                      N                      -                      1                                        N                                    ·                  N                  ·                  A                                +                A                            =                                                                                          N                      -                      1                                        N                                    ·                  B                                +                A                                                                        (        14        )            
The above equations can be used to derive another formula for the recursive imaging:
Sl(n)(t)=c1xc2x7Sl(nxe2x88x921)+c2xc2x7slk(n)(t)xe2x80x83xe2x80x83(15)
where c1 and c2 are weighting coefficients. The constants must be chosen to ensure a constant signal level and the following condition must be fulfilled: Sl(n)(t)=Sl(nxe2x88x921)(t).
The difference between equations (12) and (15), is that instead of being subtracted the information obtained by emitting with element k decays exponentially with time. In this way the information from the past is less prone to introduce motion artifacts in the image. The other benefit is that less memory is needed, since only two frames are stored.
Equation (15) can be rewritten as a sum of the signals from two consecutive emissions:
Sl(n)(t)=xcexa3bnxe2x88x92mxc2x7slk(m)(m)(t)+slk(n)(n)(t)xe2x80x83xe2x80x83(16)
assuming that c2=1 and c1=bxe2x89xa61. Nxmt less than n less than ∞ is the number of the frame, which is equal to the number of the current emission. The same element is used after every Nxmt emissions. Equation (16) can be expressed as a sum of the contributions from the emissions with the different elements:                                           S            l                          (              n              )                                ⁢                      xe2x80x83                    ⁢                      (            t            )                          =                              ∑                          k              =              1                                      N              xmt                                ⁢                      xe2x80x83                    ⁢                                    C              lk                              (                n                )                                      ⁢                          xe2x80x83                        ⁢                          (              t              )                                                          (        17        )            
Here contribution of one element is the amount of information Clk(n), added to the signal Sl(n) by the emissions with element k up to the moment of emission n. The contribution of the current emitting element (k(n)=i) to Sl(n) is:
Cli(n)(t)=sli(n)(t)+bNxmtxc2x7sli(nxe2x88x92Nxmt)(t)+b2Nxmtxc2x7sli(nxe2x88x922Nxmt)(t)+xe2x80x83xe2x80x83(18)
This is a geometrical progression. If the tissue is motionless then:
sli(n)(t)=sli(nxe2x88x92Nxmt)(t)= . . . =sli(t)
                                          C            li                          (              n              )                                ⁢                      xe2x80x83                    ⁢                      (            t            )                          =                                            [                              1                +                                  b                                      N                    xmt                                                  +                                  b                                      2                    ⁢                                          N                      xmt                                                                      +                …                            ]                        ⁢                          xe2x80x83                        ⁢                          s              li                        ⁢                          xe2x80x83                        ⁢                          (              t              )                                =                                                    s                li                            ⁢                              xe2x80x83                            ⁢                                                (                  t                  )                                ·                                                      ∑                                          p                      =                      0                                        ∞                                    ⁢                                      xe2x80x83                                    ⁢                                      b                                          pN                      xmt                                                                                            =                                          s                li                            ⁢                              xe2x80x83                            ⁢                                                (                  t                  )                                ·                                  1                                      1                    -                                          b                                              N                        xmt                                                                                                                                                    (        19        )            
If b=0.8 and Nxmt=64 then Cli(n)(t)≈sli(n)(t). This means that at emission n the contribution from the emitting element is approximately equal to the newly acquired information.
With the chosen values of b and Nxmt the contribution to the signal Sl from the previous emission with the same element is 128 dB lower than the contribution from the current one and motion artifacts can be neglected.
Improving the Signal-to-noise Ratio
The advantage of the above presented imaging approach is that a dynamically focused image in transmit and receive is obtained. This is possible because a single small transducer element emits an almost spherical wave, which unambiguously determines the propagation time of the ultrasound energy.
The drawback is that the energy, which can be sent into the body with only one array element is not enough for obtaining a high penetration depth. One way to increase the penetration depth is to emit energy with several elements as suggested in [2]. The delays of the elements are set in order to approximate the radiation pattern of a single element. Here the delays are calculated by:                               d                      i            a                          =                                            R              xmt                        -                                                            R                  xmt                  2                                -                                                      (                                                                  x                                                  i                          a                                                                    -                                              x                        a                                                              )                                    2                                                              c                                    (        20        )            
where a is the index of the central element, xa is the coordinate of the central element of the active aperture, xia is the coordinate of the element whose delay in transmit dia is calculated, and Rxmt is the radius of the spherical wave. If the number of elements in the active aperture is Nactive, then Rxmt is selected to fulfill the condition:
Rxmtxe2x89xa7(Nactivexe2x88x921)/2xc2x7pitchxe2x80x83xe2x80x83(21)
where pitch is the distance between the centers of two neighbor elements in the array. The elements"" indices are relative to the index of the central element and are given by:
ia=xe2x88x92(Nactivexe2x88x921)/2 . . . (Nactivexe2x88x921)/2xe2x80x83xe2x80x83(22)
For convenience Nactive is usually an odd number to ensure the presence of a central element. The delays in receive are calculated by formula (8), assuming that the transmitting element has the coordinate xa.
The result of using multiple elements in transmit is that the signal-to-noise-ratio is increased. Let the number of transmissions be Nxmt and the number of receiver elements Nrcv. Then the signal-to-noise ratio is proportional to:
SNRxcx9c{square root over (Nxmtxc2x7Nrev)}xe2x80x83xe2x80x83(23)
If the number of the elements in the active sub-aperture is Nactive then the SNR becomes proportional to:
SNRxcx9c{square root over (Nativexc2x7Nxmtxc2x7Nrev)}
Thus, using an active aperture in transmit with Nactive=11 results in an 11 dB increased in SNR.
By using this approach a continuous image consisting of a number of simultaneously beam formed lines can be made at the pulse repetition frequency of the emitted ultrasound, which makes it possible to follow tissue motion.
By using a array transducer a sector scan image can be made.
By using a matrix transducer a volumetric image can be made.
By using several transducer elements during the transmission an improved signal-to-noise ratio can be obtained.
By using recursive imaging velocity images can be improved and velocity distributions can be found throughout the image.